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ABSTRACT 

Detecting redshifted 21cm emission from neutral hydrogen in the early Universe promises 
to give direct constraints on the epoch of reionization (EoR). It will, though, be very chal- 
lenging to extract the cosmological signal (CS) from foregrounds and noise which are orders 
of magnitude larger. Fortunately, the signal has some characteristics which differentiate it 
from the foregrounds and noise, and we suggest that using the correct statistics may tease out 
signatures of reionization. We generate mock datacubes simulating the output of the Low Fre- 
quency Array (LOFAR) EoR experiment. These cubes combine realistic models for Galactic 
and extragalactic foregrounds and the noise with three different simulations of the CS. We fit 
out the foregrounds, which are smooth in the frequency direction, to produce residual images 
in each frequency band. We denoise these images and study the skewness of the one-point dis- 
tribution in the images as a function of frequency. We find that, under sufficiently optimistic 
assumptions, we can recover the main features of the redshift evolution of the skewness in the 
21cm signal. We argue that some of these features - such as a dip at the onset of reionization, 
followed by a rise towards its later stages - may be generic, and give us a promising route to 
a statistical detection of reionization. 
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1 INTRODUCTION 

Between redshift z ~ 20 and z « 6 the Universe underwent a 
trans ition from being almost entirely neutral to almost en tirely ion- 
ized teenson et aLil2006l: Ipurlanetto, Oh & Briggsll2006h . This pe- 
riod, the epoch of reionization (EoR), saw the first collapsed objects 
emit radiation which heated and ionized the surrounding diffuse 
gas (additional sources of heating and ionization, such as dark mat- 
ter decay, have also been considered; see, e.g.. ,Valdes et al.,,2007i) . 
Studying the emission from thi s gas, in particular the redshifted 
21cm line of n eutral hydrogen ('FielJIigSS', I1959I; 'Hogan & Ree3 
19791; 



Scott & Rees 1990; Kumar, Subramanian & Padmanabha 



I995I ; iMadau. Meiksin & Reesll 19971) . may therefore tell us about 
the physics of these objects, and about structure formation in a red- 
shift range that has previously been explored rather less directly. 
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For example, quasar absorption spectra can constrain the proper- 
ties of the inter galactic medium towards the end of reoinization 
jFan et alj2006h . while measurements of temperature and polariza- 
tion anisotropics in the cosmic microwave background provide an 
integral constraint on the density of free ele ctrons between the o b- 
server and the surface of last scattering (e.g. Dunklev et al. 2008). 

Several current and upcoming facilities (e.g. GMRtQmWaQ 
LOFArQ 21CMa|3 PAPErQ SKa|!1) win be sensitive to emis- 
sion of the right wavelength to detect a signal from neutral hydro- 
gen during the EoR. Significant observational challenges must be 
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overcome, however, before a convincing detection can be made. 
The Galactic and extragalactic foregrounds have a mean amplitude 
around 4-5 orders of magnitude larger than the expected EoR sig- 
nal (though their fluctuations, which are the relevant quantity for an 
interfero meter, are only arou nd three orders of magnitude larger; 
see, e.g., IShaveretal.lll999t) . Even closer to home, the signal is 
corrupted by the ionosphere, radio frequency interference and in- 
strumental effects. Assuming all these factors can be dealt with, for 
realistic integration times with the imminent generation of facilities 
the random noise on the measurement per resolution element will 
still be a few times larger than the signal. 

The prospect of new observations of a poorly con- 
strained period in the Universe's history also poses a chal- 
lenge to theorists: to model and characterize the 21cm emis- 



sion in such a way that it can be meaningfully com- 
pared to the data (e.j' ^ — ' ' — ^ 'onm'. t r> — i 



Barkana & Loeb 2001; Loeb & Barkana 
2001l;ICiardi. Ferrara &. White 2003; Ciardi, Stoehr & White 2003; 



Bromm & Larson 2004; Gleser et al. 2006; Ilie v et al.ll200A 120081 : 
Zaroubi et al..,2007. ; , Thomas & Zaroubi,2008) . Given the observa- 



tional limitations listed above, it is unlikely that there will soon 
be clean maps of the EoR signal with which to confront mod- 
els. The first detection of reionization from redshifted 21cm data 
will therefore be of a statistical nature. This raises the question 
of precisely which statistics to use: a questi on discussed by, e.g., 
iFurlanetto. Zaldarriaga & Hemguistl (l2004al lbh; iBharadwai & Alii 
l l2005h ; iGleser. Nuss er & BensonI ( l2008h . In the first instance, they 
should provide a clear indication of the global transition from a 
Universe that is mostly neutral to one that is mostly ionized. Ide- 
ally, they should be able to discriminate between different models 
describing the more detailed progress of reionization. Most impor- 
tantly, though, they should be robust to the contamination intro- 
duced by the observing process, and to the presence of high levels 
of noise. 

We propose using the skewness of the one-point distribution 
of the brightness temperature to study reionization, though we also 
consider the prospects of other, similar statistics: the unnormalized 
third moment and the kurtosis of the distribution. As we shall see 
below, general arguments suggest that the skewness should be a 
strongly evolving function of redshift during the EoR, and these 
arguments are supported by simulations. Using these simulations, 
we generate datacubes which also incorporate realistic models for 
the foregrounds, instrumental response and noise levels expected 
for the LOFAR EoR experiment. We generate residual images at 
each observed frequency by attempting to remove the foregrounds 
using a fitting algorithm, then study the properties of these residual 
images as a function of redshift. If these residual images are de- 
noised appropriately, we find that we can indeed track the progress 
of reionization using the skewness. 

In Section |2] we introduce the skewness and explain how it 
may help. In Section[3]we give a brief description of our models for 
the cosmological signal (CS), instrumental response, foregrounds 
and noise. Then, in Section|4l we describe our method for extract- 
ing the signal from datacubes which combine all these components, 
and present our results. We discuss some possible problems with 
and extensions to our methods in Section |5] and finally we offer 
some conclusions in Section|6] 



2 GENERAL APPROACH 

We assume, with reasonable observational support, that the fore- 
grounds are smooth as a function of frequency, and exploit this in 



order to extract th e CS dShaver e t al.ll 199 91; loi Matteo et al.ll2002l ; 
lOh & Mackl |2003| ; IZaldarriaga. F urlanett o & Hemguistl l2004h . A 
simple way to imagine doing this is first to estimate the fore- 
grounds, either by fitting a smooth function to the intensity as a 
function of frequency along each line of sight, or by applying a 
filtering procedure. Subtracting this estimate from the total yields 
fitting residuals which are an estimate of the CS plus random noise. 

If the variance of the noise is known, then at each frequency 
(or equivalently at each redshift) we can attribute excess variance 
in the fitting residuals over and above this level as coming from 
the CS. This may yiel d a statistic a l dete ction of emission from the 
epoch of reionization. Ijelic et al] j2008l) demonstrate this method 
using the same foregrounds and noise levels as are used in this 
work. 

There are some difficulties with both parts of the above pro- 
cedure. Firstly, we may over-fit or under-fit, leading (respectively) 
to underestimating or overestimating the residuals. This problem is 
exacerbated near the ends of the frequency range. Secondly, the 
noise properties may not be sufficiently well characterized. The 
variance of the signal is expected to be only a small fraction of 
the variance of the noise, and hence the latter must be very well 
known. 

We therefore seek a statistic on the fitting residuals which will 
detect the onset of reionization more robustly with respect to er- 
rors in our estimates of the variance. A possible candidate is the 
skewness, 71, defined in general for a continuous distribution with 
probability density function f(x) as 



71 



^3 



(1) 



where a is the variance of the distribution and /is is the third mo- 
ment about its mean, ^. A distribution which is mainly concentrated 
at low X, but with a tail towards high x, will have positive skew- 
ness. Similarly, a distribution with a tail extending to low x will 
have negative skewness, and a poor estimate of a cannot change its 
sign. 

Rather than dealing with a continuous distribution, we will be 
computing the skewness for images with A'^ pixels, in which the ith 
pixel has a temperature Ti. Then the skewness may be expressed as 



71 = 



(2) 



where T is the mean temperature in that image and the sums are 
over all pixels. In the case of residual images, if the foreground 
fitting is unbiased, and if the noise is not skewed, then any signifi- 
cant skewness remaining must come from the cosmological signal. 
Moreover, we may expect some skewness in the cosmological sig- 
nal, which becomes very non-Gaussian once ionized bubbles ap- 
pear in large numbers during reionization. 

To be concrete, making some reasonable assumptions and ap- 
proximations (that the optical depth is much less than unity, that the 
spin temperature of the neutral hydrogen is much greater than the 
CMB temperature, and that at these redshifts the Hubble parame- 
ter H{z) « n]l^HQ{l + zf'\ the difference, 5T^, between the 
bri ghtness t emperature of the 21c m emission and the CMB is given 
bv JMadau et al. 1997i; ICiardi & Madau 2003) 



57b 
mK 



= 39/1(1 + (5)XHI 



0.042 



0.24 



10 



(3) 



where 5 is the matter density contrast, xm is the neutral hy- 
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Figure 1. The distribution of 5T^ in the f250C simulation (see Section [3Tl 
at four different redshifts, showing how the distribution evolves as reioniza- 
tion proceeds. Note that the y-axis scale in the top two panels is different 
from that in the bottom two panels. The delta-function at (5Tb = grows 
throughout this period while the rest of the distribution retains a similar 
shape. The bar for the first bin in the bottom-right panel has been cut off: 
approximately 58 per cent of points are in the first bin at z = 7.78. 



drogen fraction, and the current value of the Hubble parameter, 
Ho = lOO/i kms"^ Mpc-\ 

At high redshift, when a::Hi is close to unity everywhere, the 
distribution of intensities is governed by the density field, 1 + 5. Ini- 
tially this is nearly Gaussian, but develops a positive skewne ss due 
to gravitational instabilities: see, for example, IPeebiej l ll980l) . This 
period is illustrated in the top left panel of Fig.[T] which shows the 
one-point distribution of STh in one of our simulations (the f250C 
simulation; see Section lJTt at z — 10.6, corresponding to an ob- 
served frequency of 122.5 MHz. If reionization then takes place in 
patches, with large volumes remaining mostly neutral while almost 
fully ionized bubbles form around sources of ionizing photons, this 
has the effect of setting xm ~ (and so ST^ — 0) within the 
bubbles, affecting the distribution of STh outside the bubbles only 
weakly. So, in an idealized case, reionization takes points from the 
distribution of ST^, and moves them to a Dirac delta-function at 
zero. This has the effect of making the skewness less positive; it 
may even become negative. The distribution of ST^ at an early 
stage in this process (z — 9.12, or 140.3 MHz) is shown in the 
top right panel of Fig.Q By z = 7.98 (158.2 MHz; bottom-left 
panel) the two parts of the distribution are very distinct. This may 
help to make it clear how the skewness can vanish: the mean ST^, 
lies between these two peaks, and the negative contribution to the 
skewness from points in the delta function at zero may cancel with 
the positive contribution from points to the right. At a later stage 
of reionization, when most of the pixels in a noiseless map of STh 
at a given frequency have values near zero, the points outside ion- 
ized bubbles form a high-^Tb tail, giving the overall distribution a 
strong positive skew. This can be seen in the bottom right panel 
of Fig. 0(2 = 7.78 or 161.8 MHz). Note the short time be- 
tween the third and fourth panels: the later stages of reionization 
can progress rather quickly as the number of ionizing sources can 



rise very rapidly, especially if a major part is played by massive 
sour ces residing in haloe s in the exponential tail of the mass func- 
tion jjenkins et alj200lh . 

In this idealized case, then, the skewness as a function of red- 
shift should show a dip in the early stages of reionization, before 
growing large in the later stages. Our aim in the subsequent sec- 
tions of this paper is to test if such a characteristic feature is indeed 
seen in realistic simulations of reionization, and whether or not it 
can provide a robust detection; or in other words, whether the ef- 
fects of foreground subtraction, noise, and instrumental corruption 
can mask or mimic the signal. 



3 SIMULATIONS 

3.1 Cosmological signal 

We use three simulations to estimate the CS. The first a n d mos t 
detailed is the simulation labelled f250C by llliev et al.l ( l2008h . 
The methodology behind t his simulation is mor e fully described 
by Iliev et al. ( 2006) and Mellema et al. (2006lj |). The cosm olog- 
ical particle-mesh code PMFAST (Merz. Pen & Trad 12005) was 
used to follow the distribution of dark matter, using 1624^ par- 
ticles on a 3248^ mesh. The ionization fraction was then cal- 
culated in post-processing using the radiative transfer and non- 
equilibrium chemistry code C^-RAY iMellema et al.ll2006ij) . This 
takes place on a coarser, 203'^ mesh, and this is therefore 
the size used in this work. The simulation box has a comov- 
ing size of 100 Mpc. The cosmological parameters are 

close to those inferred from the thr ee-year Wilkinson M icrowave 
Anisotropy Probe data (WMAPy. ISpergel et aEI l2007h . namely 
(fim, ^A, ^b, h, as, n)=(0.24, 0.76, 0.042, 0.73, 0. 74, 0.95). 

A slightly different approach, detailed by iThomas et al.l 

( l2009l) . is used to generate our other simulations. The dark 
matter distribution is calculated using the TREE-P M A^-body 
code GADGET2 dSpringel. Yoshida & White! 1200 ll; ISpringell 
l2005l) . Ionization is the n calculated using a on e-dimensional 
radiative transfer code jThomas & Zaroubil |2008|) . The speed 
of this approach means it is possible to study many more al- 
ternative scenarios for the reionization process, while retaining 
good agreement with more accurate, three-dimensional radiative 
transfer simulations. We will show results from two different 
simulations. In both cases, the dark matter simulation uses 512'' 
dark matter particles in a box of comoving size 100 /i"^ Mpc, with 
{Qn,, f^A, f^b, h, (78, n)=(0.238, 0.762, 0.0418, 0.73, 0.74, 0.951). 
While the simulations contain no baryons, this value of f2b was 
used to generate the initial power spectrum. These parameters 
give them lower resolution than the dark matter part of the f250C 
simulation, meaning that low mass sources are not resolved and 
are neglected. In one of these simulations we assume that the 
Universe is reionized by QSOs, and in the other by stars. These 
two simulations are labelled 'T-QSO' and 'T-star' respectively. 
The former should not be affected too seriously by the lack of 
resolution, since QSOs do not reside in low-mass haloes. In the 
latter, the geometry of reionization may be altered: compared to 
a higher resolution simulation, larger ionized bubbles may form 
at a given global star formation rate, for example. As we shall 
see below, 'T-star' shows rather different characteristics from the 
f250C simulation, despite the fact that stars provide the ionizing 
photons in both cases. This illustrates the uncertainties involved 
in modelling the physics of reionization, in selecting the source 
populations and finding their distribution in space, and in choosing 
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the approximations required to make the calculations tractable. We 
do not analyze the differences between the simulations in great 
detail here; rather, we use the different simulations to provide 
a variety of plausible scenarios with which to test our signal 
extraction techniques. 

The above calculations all take place in periodic boxes. The 
final step in generating a datacube is to take a series of simula- 
tion snapshots at different redshifts, and interpolate between them 
to produce a spatial slice at each obse rved frequency . This pro- 
cedure, which is desc ribed in detail by iThomas et al.l ( l2009h and 
lMellemaetd] ( l2006b! ), is analogous t o the gene ration of lightcone 
output to compare to galaxy surveys jEvrard e t al. 2002). In per- 
forming this conversion from a position in a periodic box to a red- 
shift, we take account of the peculiar velocities; that is, our dat- 
acubes are in redshift sp ace, as will be the case f or the observational 
data. As emphasized bv lMellema et alj j2006bh . who were the first 
to include peculiar velocity distortions in a redshifted 21cm con- 
text based on detailed simulations, these effects can be important. 
At linear scales, the redshift space distortions have th e effect of en- 
hancing density fluctua tions along the line of sight ( iKaiserll 19871 ; 
iBharadwai & Aiill2005h . 

3.2 Foregrounds, noise and instrumental effects 

We us e the foreground simulations described in detail bv ljelic et al.l 
( l2008h . These incorporate contributions from Galactic diffuse syn- 
chrotron and free-free emission, and supernova remnants. They 
also include extragalactic foregrounds from radio galaxies and ra- 
dio clusters. The foreground maps cover an area 5° x 5° on the 
sky, which corresponds to the area of one LOFAR EoR w indow. We 
also ad opt the frequency-dependent noise levels given bv lJelic et al.l 
( l2008h . The noise is described in more detail below, when we in- 
troduce each of our two noise models. 

In Section |4] we will give results obtained by combining the 
CS, foregrounds and noise, and then attempting to extract the CS. 
We also attempt the more difficult, and more realistic task of ex- 
tracting the signal given such a datacube corrupted by the instru- 
mental response (so-called 'dirty maps '). These instrumental ef - 
fects will be described in more detail bv lLabropoulos et aLll l2009l) . 
In brief, at each observed frequency, each pair of LOFAR stations 
gives an estimate of the Fourier transform of the sky brightness 
along a track of points (a uv track) in Fourier space (the uv plane). 
Translated into configuration space, this means that our images of 
the sky are convolved with a complicated point spread function (the 
'dirty beam'). Moreover, because the origin of the Fourier plane is 
not sampled, interferometer measurements are not sensitive to the 
mean brightness. A more detailed introduction to such datacubes, 
and their r elation to the power spectr um of emission from the EoR, 
IS given bv lMorales & HewittI ( [20041) . 



4 RESULTS 

4.1 Skewness in perfect data 

The evolution of skewness in our three simulations, uncorrupted by 
noise or foregrounds, for the redshift range observable by LOFAR, 
is given in Fig. [2] The most obvious feature here, common to all 
three simulations, is the rise in the skewness at low redshift, due to 
the high-^Tb tail of points with some remaining neutral hydrogen. 
At higher redshift there are some differences, however. In f250C 
there is a clear dip &X z ^ 7.8-9, whereas this is not so obvious in 




Figure 2. The evolution of the skewne.ss of the distribution of 5T^, as a 
function of redshift in our three simulations. Each point corresponds to the 
skewness of the one-point distribution in a slice through a datacube at the 
given redshift. The most striking feature here is the steep rise in skewness at 
low redshift. The details at high redshift can be seen more clearly in Fig. [3] 



either of the other simulations. For the T-QSO simulation, though, 
and more noticeably for the T-star simulation, the skewness is nega- 
tive in certain redshift slices. All the simulations show large spikes 
in skewness at low redshift. At these redshifts, a slice may have 
only one or a few regions from which there is significant emission. 
This leads to large variation between slices. A slice in which one 
small region produces a high-emission tail in the one-point distri- 
bution, and which is surrounded by more uniform slices, shows up 
as a spike in the evolution of the skewness. 

The lower-skewness region of the plot is shown in more de- 
tail in Fig. [3] where we also compare the evolution of the skewness 
to that of the mean. Whilst the mean value of 5T\^ is unobserv- 
able with an interferometer such as LOFAR, this serves to illus- 
trate the progress of reionization. Clearly the process is much more 
rapid in f250C than in the T-QSO simulation, which in turn is more 
rapid than the T-star simulation. In Fig. [3] a dip in the skewness 
is discernible in T-QSO, but is much less obvious than the dip in 
f250C because of fluctuations at high redshift, and because it spans 
a wider range in redshift due to the more extended reionization. 
This may be an indication that a more extended reionization process 
will be harder to detect using the skewness. The T-star simulation 
provides some hope, however. Despite an even more gradual reduc- 
tion in the mean differential brightness temperature, the skewness 
is negative for quite a large range in redshift. Because the density 
field is positively skewed, a negative skewness is a clear signature 
of reionization. It is therefore possible that the skewness could pro- 
vide a detection even in the case of very extended reionization. 

4.2 Extracted skewness 

We now proceed to test the possibilities for signal extraction using 
the skewness, starting with a rather more optimistic case than will 
be encountered with the actual LOFAR EoR experiment. We first 
note that the 5° x 5° field of one LOFAR EoR window corresponds 
to a distance of approximately 800 Mpc x 800 Mpc (comoving) 
at 2 = 10 in our assumed cosmology. At each redshift we there- 
fore tile copies of the simulation to produce a slice of the correct 
size, then interpolate this onto a 256^ mesh. Since each pixel will 
be affected differently by foregrounds and noise, and since we con- 
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Figure 3. The evolution of the mean differential brightness temperature (top 
panel) and of the skewness (bottom panel) in the three simulations. The 
colours and styles of the lines are the same as in Fig. [2] The top panel 
shows that reionization is extended to very different degrees in the different 
simulations. We cut off the large values of the skewness in the bottom panel 
so that the dip in skewness at high redshift can be seen more clearly than in 

Fig.m 



sider only one-point statistics, we do not anticipate that this will 
strongly affect our conclusions. To produce our datacube, we add 
the simulated foregrounds described in Section [T2l to the CS, then 
smooth each slice using a Gaussian kernel to the estimated LO- 
FAR resolution of ^ 4 arcmin. We the n add uncorrelated random 
noise as described bv lJelic etalj j2008l) . The noise has an rms of 
52 mK at 150 MHz and has two contributions: a frequency de- 
pendent component coming from the sky, which scales as u^'^'^^, 
and a frequency-independent part from the receivers. The noise on 
each image pixel is independent. In reality, this will not be the case: 
rather, the noise on individual visibilities will be independent. We 
will tackle this more difficult case with realistic noise and a non- 
Gaussian point spread function below. Spatial slices are separated 
by 0.5 MHz in frequency, v. At 150 MHz this corresponds to a 
difference in redshift, Az ~ 0.03, the slices having a comoving 
thickness of approximately 7 Mpc. 

Once we have a datacube with EoR signal, foregrounds and 
noise, we fit a third-order polynomial in log u to each pixel. We 
have experimented with using different functional forms, but find 
that so long as we obtain a visually reasonable fit, our results for 
the skewness do not change enough to affect our conclusions. Mea- 
surements of the variance are rather more sensitive to under- and 
over-fitting, which demonstrates the importance of understanding 
the foregrounds well, and of using robust statistics. It is also possi- 
ble to estimate the foregrounds by removing noise using a filtering 
procedure. While this requires fewer assumptions about the nature 
of the foregrounds, it tends to over-fit. 

After a fit has been obtained, it is subtracted from the total, 
leaving residuals which are an estimate of the CS plus the noise. 
The skewness of these residuals as a function of redshift is shown 
in Fig.|4l While the CS which goes into the datacube is different 
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Figure 4. Skewness of the fitting residuals from datacubes with uncon'e- 
lated noise. The same noise is used in each cube. The colours and styles 
of the lines are the same as in Fig. [2] in the top panel, the component of 
the datacube from the CS comes from f250C; in the middle panel, from the 
T-QSO simulation; and in the bottom panel, from the T-star simulation. The 
similarity between the three panels arises because the noise realization (and 
the foregrounds) is the same for each panel, and dominates over the CS. 



for each panel of the plot, the noise and foregrounds are the same. 
This accounts for the fact that the skewness as a function of redshift 
shows very similar features in each panel. At first sight this seems 
rather discouraging, with the desired signal totally dominated by 
fitting errors (FEs) and noise. 

The situation can be improved, however. At each frequency, 
we may denoise the residual image by smoothing. This is possible 
because our images are oversampled, with more than one pixel per 
resolution element, and because there are no pixel-to-pixel corre- 
lations in the noise (by construction). While this is clearly unreal- 
istic, it serves as a prototype for the more difficult denoising step 
when we consider the dirty maps. The effect of smoothing on the 
different components of the residual maps - these components be- 
ing the CS, FEs and noise - is illustrated in Fig. [5] in which we 
show how the absolute value of the third moment of the one-point 
distribution of these components, \fi3\, is affected when they are 
smoothed with windows of different size. We show the result for 
the slice of the T-star datacube at 115 MHz, corresponding to a 
redshift of 11.35. This slice is chosen because the skewness of the 
T-star simulation is significantly negative here; we get similar re- 
sults with the other simulations if we choose an appropriate slice 
in which the skewness is significantly different from zero. When 
the smoothing window is very narrow, so that there is almost no 
smoothing, Ifial for the noise exceeds the value for the CS. This 
occurs even though (/i""""^) — (where the expectation is taken 
over different noise realizations), simply because the noise rms is 
so much larger than that of the (significantly skewed) cosmological 
signal. As the size of the smoothing window is increased, Iais"'""! 
drops much more quickly than l/ig ^ | since the smoothing averages 
together uncorrelated noise pixels, but correlated signal pixels. At 
large scales, the signal also becomes uncorrelated, so its small rms 



means that \^^^ \ 



I < 1^3°'"°! once more. The scale at which 
by the greatest amount in this residual map is ap- 
proximately 3-4 arcmin. In the case of the fitting errors, \fi3^\ 
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Figure 5. The effect of smoothing scale on the absolute value of the third 
moment of the one-point distribution, \ (defined in Equation[T), for the 
three different components of the residual maps at 115 MHz (z = 11.35). 
The width, dg, of the Gaussian distribution which forms our smoothing 
kernel is defined here as the distance between x = iter for a distribution 
with standard deviation cr. The CS (blue, solid line) is from the T-star sim- 
ulation, in which the skewness of the one-point distribution is negative at 
this frequency. The red, dashed line corresponds to the noise (where here 
the noise on each pixel in the unsmoothed map is independent), and the 
green, dot-dashed line to the fitting errors, by which we mean the difference 
between the foregrounds cube and the polynomial fit to the full datacube. 
The value for the foregrounds themselves, before fitting, is several orders of 
magnitude larger 
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Figure 6. Skewness of the fitting residuals from datacubes with uncon'e- 
lated noise, but in which the residual image has been denoised by smooth- 
ing at each frequency before calculating the skewness. The three lines cor- 
respond to the thr'ee simulations, with colours and Hne styles as indicated in 
Fig. [2] Each line has been smoothed with a moving average (boxcar) filter 
with a span of nine points. The grey, shaded area shows the errors, esti- 
mated using 100 realizations of the noise. The fact that the lines differ so 
much is in marked contrast to Fig. |4] and shows the impact of smoothing 
the residual images to suppress the noise. 



shows less variation as the smoothing scale is changed than either 
of the other components. If the foreground fitting were completely 
unbiased, we might expect that any errors in the fitting would be 
Gaussian and caused entirely by noise, and so this component of 
the residual images would behave similarly to the noise compo- 
nent. The fact that the skewness in the fitting errors appears to come 
partly from large scales suggests that bias in the fitting may allow 
some leakage through from the foregrounds themselves, which are 
correlated on large scales. If the skewness of the foregrounds is 
larger than we have assumed, therefore, we will need to fit them 
more accurately in addition to exploiting this scale dependence. In 
the present case, j/i™] is similar to |^3°"^°| at the scale at which 
the latter is dwarfed by the contribution from the CS. 

In practice, to extract the skewness as a function of frequency 
we make the natural choice of smoothing scale, using the same ker- 
nel as was used to degrade the signal and foregrounds to the resolu- 
tion of the telescope. We then compute the skewness in each slice as 
before. The skewness as a function of frequency for each datacube, 
after following this procedure, is shown in Fig. |6] To improve the 
clarity of the plot, each line is smoothed by taking a moving av- 
erage with a span of nine points (a boxcar filter). To estimate the 
error, we generate 100 datacubes containing the foregrounds and 
with different realizations of the noise, but with no CS present. We 
feed each cube through our fitting and smoothing procedure, calcu- 
late the skewness as a function of redshift, and smooth this function 
with a moving average filter just as for the cubes containing a sig- 
nal. The range between the 16th and 84th percentile of the skewness 
for these realizations is shown as the light grey shaded area in the 
figure. 

One can see from Fig.[6]that this smoothing procedure allows 
us to extract a significant signal, despite making only rather general 



assumptions about the scale at which features due to the signal, in- 
strument and noise are important. The result for f250C (blue, solid 
line) is most striking, with rapid transitions in the skewness in the 
range z ~ 7.5-8.5. The position of the dip corresponds to the posi- 
tion of the dip in the uncoirupted simulation shown in Fig. [2] While 
the skewness continues to rise in the original simulation, however, 
for the extracted signal it returns to zero at low redshift. This is be- 
cause the variance of the CS becomes very small at low redshift. In 
the uncorrupted simulation, this allows the skewness to grow very 
large. In the residual images, however, the variance of the noise and 
fitting residuals comes to dominate, even after smoothing, which 
drives the skewness towards zero. We return to this point below 
when we consider alternative statistics. The extracted signal for the 
other two simulations shows the behaviour one might expect: the T- 
QSO simulation (red, dashed line) shows only a weak dip in skew- 
ness, but a strong peak due to the rapid rise in skewness for the 
uncorrupted simulation at 2: < 8.5. The T-star simulation, mean- 
while, shows a gradual rise in skewness throughout the redshift 
range, with significant non-zero skewness detected for 2 > 9.5 
and z < 7.5. 



4.3 Skewness from dirty maps 

We now move on to an analysis of the so-called 'dirty maps'. To 
generate these, we first add together the unsmoothed foreground 
and signal cubes, where the latter have been tiled, as before, to pro- 
duce maps of the right angular size. Each slice is then corrupted by 
the instrumental response. We achieve this in practice by Fourier 
transforming the image, multiplying by the sampling function (cal- 
culated on a grid with the same number of points as the image), 
and then applying the inverse Fourier transform. This is equivalent 
to convolving each image with the point spread function (PSF) of 
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the instrument. We make the simplifying assumption that the sam- 
pling function does not change with frequency. If the uv plane is 
uniformly filled, this should not be excessively optimistic. It could 
be enforced in practice by discarding high-fc data so that equiva- 
lent (and completely filled) areas of the uv plane are retained in 
each frequency band. The noise is dealt with slightly differently. 
We consider pixels in the uv plane where the sampling function is 
non-zero to be encompassed by our uv coverage, and we generate 
uncorrelated Gaussian noise at each such pixel. Pixels outside our 
uv coverage are set to zero. We (inverse) Fourier transform to re- 
turn to the image plane, then normalize this 'noise image' such that 
it has the correct rms. This procedure yields noise that is almost 
uncorrelated between independent resolution elements (though the 
noise on adjacent pixels is correlated). 

We fit out the foregrounds in the dirty cubes in the same way 
as before. The skewness of the residual images exhibits the same 
problem seen in Fig.|4l being dominated by the noise. In this case, 
the smoothing procedure used above would not be expected to help, 
since the noise is correlated on the scale of our smoothing kernel. 
In addition, since our resolution is comparable to the scale of fea- 
tures in the original signal, using a broader kernel simply washes 
out the signal as well as the noise. We therefore require a more 
sophisticated denoising scheme. 

4. 3. 1 Wiener filtering 

With the results of Section l42l in mind, we use the differing cor- 
relation properties of the signal and noise in our extraction. To be 
explicit, suppose that we write the residuals as a vector d, where 
di is the residual at the ith pixel of a map at a given frequency. We 
relate d to the image from the uncorrupted simulation, s, by 



d = Rs + e 



(4) 



The matrix R encodes the convolution of the signal with the PSF, 
while e represents the noise. We neglect any contribution to e com- 
ing from errors in the fitting procedure, so we can assume that the 
correlation matrix of the noise, N = (ee^), is known (here, is 
the conjugate transpose of e). 

We consider a very optimistic situation for extracting the 
skewness, which occurs if the correlation matrix of the signal, 
S = (ss^), is also known. We can then perform a Wiener decon- 
volution on each residual image to recover an estimate of the CS. 
That is, we compute s = Fd where the Wiener filter F is given by 

F = SRt(RSRt (5) 

(see, e.g., IZaroubietal.lll995l) . In the absence of noise, this pro- 
cedure reduces to an ideal inverse filter that estimates the orignal 
image before corruption by the PSF. In the presence of noise, the 
Wiener filter suppresses power in the image at those values of k for 
which the signal-to-noise ratio (SNR) is low , while retaining p ower 
for modes where the SNR is high (see, e.g.. |Press et al.lll986l) . The 
algorithm is optimal in the least-squares sense. 

The skewness of these deconvolved images as a function of 
redshift is shown in Fig.|7] Comparing to Figs.|2]and[3] one can see 
that this procedure gives excellent results, recovering the general 
trends in skewness seen in the original simulations. Indeed, using 
an optimal filter with precise knowledge of the signal and noise 
properties means that we recover larger values for the skewness 
than were seen after applying the simple smoothing to the uncorre- 
lated noise case of Section |42l (Fig.[6l). We can realistically expect 
a situation intermediate between the results of Figs. |4] and [T] The 
lines representing f250C and the T-QSO simulation do not extend 
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Figure 7. Skewne.ss of deconvolved images as a function of redshift. The 
dirty maps are generated by convolving the CS and foregrounds with the 
instrumental response, and adding noise with a realistic coiTelation ma- 
trix. The CS is then estimated by titting out the foregrounds and applying a 
Wiener deconvolution to the residuals, assuming that the PSF and the cor- 
relation matrices of signal and noise are known. Line colours and styles are 
the same as in Fig.|2] The deconvolution (as opposed to the simple smooth- 
ing which was sufficient for the uncorrelated noise used in Fig. [6) is unstable 
when the cosmological signal becomes very small, which prevents us from 
estimating the errors in the same way as for Fig. [6] 



all the way to z = 6 in Fig. [7] At these redshifts, the variance in the 
CS is so small compared to that of the noise that the deconvolution 
becomes unstable. For the same reason, we cannot estimate the er- 
rors in the same way as for Fig.[6](that is, by generating realizations 
with no CS at all). None the less, the errors can be inferred to be 
small since we do not see the same effect as in Fig.|6l in which the 
extracted signal from the datacube generated with the f250C and 
T-QSO simulations is very similar at high and low redshift, being 
dominated by noise. 

An obvious objection to the method presented here is that if 
the correlation matrix of the CS is known, this means that we have 
already detected a signal from the EoR, so higher-order statistics 
are not required to extract it. The force of this objection depends 
on how good an estimate of the correlation matrix of the signal 
is required for the deconvolution to give an acceptable result. We 
present a test of this in Fig. [8] The three lines in the figure show the 
skewness extracted from the f250C residuals using three different 
assumption for the correlation matrix, S, used in the Wiener decon- 
volution. The solid, blue line shows, for reference, the skewness 
extracted when we use the correct correlation matrix Sf25oc(2)> 
calculated from the original simulation, in performing the decon- 
volution. The dashed, magenta line shows the extracted skewness 
when we use |Sf25oc(-2) instead. Underestimating the correlation 
matrix by a factor of two clearly has only a minimal effect on the 
extraction of the skewness. Finally the dot-dashed, cyan line shows 
the result when we use the correlation matrix of the T-star simu- 
lation. Even though the redshift evolution of the two simulations 
is very different, the dip and peak in the skewness at z ~ 7.5-8.5 
are recovered, though it is not clear that they can be easily dis- 
tinguished from the spurious variations at high redshift. This pre- 
liminary result is encouraging, but it would be preferable to use a 
correlation matrix estimated from the data themselves. For exam- 
ple, writing Equation [5] in terms of the correlation matrix of the 
residuals, D — (dd^), it becomes 
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Figure 8. The skewness as a function of redshift recovered from the f250C 
simulation by Wiener deconvolution, using three different assumptions for 
the correlation matrix, S, of the original signal. For the solid, blue line we 
use the correct correlation matrix, Sf250C (•^)- For the dashed, magenta line 
we use 5Sf25oc(^)' ™d for the dot-dashed, cyan line we use the correla- 
tion matrix from the T-star simulation. 
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have to subtract the foregrounds. Therefore it is the combination 
of bright foregrounds and structured noise which prevents us from 
extracting the skewness without some assumptions about the prop- 
erties of the signal. It will therefore be important to develop further 
our techniques for foreground subtraction, and to test how sensitive 
they are to variations in our models for the foregrounds. 



4.4 Skewness from Fourier space data 

The observable quantity for an interferometer is the radio visibility: 
a quantity that lives in Fourier space. It also seems natural, there- 
fore, to try to extract a signal from the EoR without transforming 
to the image plane first (though we note that some stages of the 
analysis which must be completed to produce the clean datacube 
from which we attempt to extract a signal, such as the subtraction 
of bright point s ources, are a lso carried out in the image plane). For 
example, Dat ta et al.l ( l2007i) have presented a formalism to search 
for bubbles in 21cm data using a statistic on the visibilities. 

Unfortunately, it is not convenient to calculate the skewness 
from Fourier space data. The case for the moments of the density 
field is well known: the variance of the overdensity field, (5^), is 
equal to the two-point correlation function evaluated at zero sep- 
aration, 5(0), which in turn is equal to an integral over the power 
spectrum of fluctuations. Similarly, the third moment {5'^) is equal 
to C(0, 0) where CC^'ij^'a) = {5{x)S{x + r\)S{x + r2)). Then 
we have 



C(o,o) 



J d'k'd'k"B{k',k",-k' -k") 



where B is the bispectrum, defined by 
(5(fci)5(fc2)5(fe3)> = 5o(ki + fca + fc3)-B(fci, fca, fca) 



(7) 



(8) 
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where 5{k) is the Fourier counterpart of S{x) and 5d is the Dirac 
delta-function. The bispectrum is a rather unwieldy object to grap- 
ple with in this context. Foreground extraction is also problematic: 
while the foregrounds remain smooth as a function of frequency at 
a given uv point, the angular scale sampled by that point is also a 
function of frequency. Further consideration of Fourier space statis- 
tics is therefore beyond the scope of this paper. 



Figure 9. Skewness of dirty maps with realistic noise, but in which no fore- 
grounds have been added. This is equivalent to a case in which we achieve 
perfect removal of the foregrounds. 



F=R^(D N)D-' . (6) 

Unfortunately this filter does not work so well in practice, since the 
convolution means that important small-scale information present 
in S is not present in D. Since estimation of the signal correlation 
matrix is intimately connected with techniques for power spectrum 
estimation, we defer further investigation to a future paper. 

4.3.2 Foreground subtraction effects 

Part of the need for sophisticated denoising techniques comes from 
the fact that imperfect fitting of the foregrounds introduces errors 
into our residual images. We illustrate this in Fig. [9l in which we 
show the skewness of the dirty maps to which we have added real- 
istic noise, but no foregrounds. The main trends in the evolution of 
the skewness are clearly visible. We have seen in Fig. |4] that even 
uncorrelated noise prevents us from recovering a signal if we first 



4.5 Alternative statistics 

The potential for using the skewness of the CS to help in the extrac- 
tion invites the question of whether other one-point statistics, such 
as the kurtosis, could also be useful. We define the kurtosis here as 
/i4 /(J* — 3 where /i4 is the fourth central moment of the distribution, 
and we subtract 3 so that a Gaussian distribution has a kurtosis of 
zero. In fact the kurtosis does evolve strongly in the signal simula- 
tions. As was the case for the skewness, it is not difficult to see why. 
While the brightness temperature traces the density field it retains a 
kurtosis similar to that of a Gaussian distribution. The formation of 
bubbles then produces a bimodal distribution with no strong central 
peak, so the kurtosis, a measure of the 'peakiness' of the distribu- 
tion, is reduced. In the final stages of reionization, the distribution 
becomes strongly peaked around zero, with a tail of points with 
strong emission, leading to a large kurtosis. The evolution of the 
kurtosis is therefore qualitatively similar to that of the skewness. 
We see all these stages in our simulations, and the evolution of kur- 
tosis as a function of redshift is shown in Fig. [TO] Unfortunately 
these trends seem to be much harder to recover than in the case 
of skewness. While we can weakly recover the dip for the case of 
uncorrelated noise in the image plane, using a similar procedure to 
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Figure 10. The evolution of the kurtosis of the distribution of ST^, as a func- 
tion of redshift in our three simulations. Some of the qualitative features are 
similar to those seen in the skewness in Fig.|2] though for slightly different 
reasons which are explained in the main text. In the definition of kurtosis 
used here, a Gaussian distribution has a kurtosis of zero. As in the lower 
panel of Fig. |3]we cut off the very large values at low redshift, in order to 
make the high redshift region clearer. 



that used for Fig. |6l we cannot recover the low-redshift peak. For 
the dirty maps, the dip in kurtosis cannot be seen even after the 
optimistic deconvolution used in Fig.|7] 

Because most of the variance in the unsmoothed residual im- 
ages comes from the noise and fitting errors, which one might ex- 
pect not to be skewed, we have investigated whether or not the un- 
normalized third moment, ^3, might perform better than the skew- 
ness, fis/a^, for signal extraction. In fact we find that they perform 
similarly, but the results for the skewness tend to be easier to inter- 
pret. This is because to calculate the errors we must still estimate 
the expected spread in from the random noise, and hence have an 
estimate of the rms of the noise. Because the noise power changes 
with frequency, the errors also become frequency dependent: the 
shaded area in a figure analogous to Fig. |6]no longer has constant 
width. 



5 DISCUSSION 

We have demonstrated that the skewness could be a useful tool for 
signal extraction in the presence of realistic overall levels of fore- 
grounds and noise. It could be, though, that some aspect we have 
not accounted for makes the process more difficult. For example, 
while the foreground maps we use have rms fluctuations of the cor- 
rect magnitude, they have rather low levels of skewness, perhaps 
unrealistically low. If, instead, the foregrounds turn out to be very 
skewed, then unless the algorithm to fit out foregrounds is unbi- 
ased, this could propagate into the fitting residuals and drown out 
the CS. In this case, however, the characteristic pattern in the skew- 
ness as a function of redshift - a dip followed by a peak - may 
allow the signal to be picked out even in the presence of residuals 
from the foregrounds. This would be exploiting the smoothness of 
the foregrounds as a fun ction of frequency once more. 

As pointed out by Ijelic et al.l UooA data constraining the 
characteristics of the foregrounds at the relevant scales and fre - 
quencies are quite scarce (though see, for example. lPen et al.l2008h . 
None the less, the extrapolations we make from larger scales and 



higher frequencies may be pessimistic, if anything. Moreover, the 
structured noise appears to be at least as influential as the fore- 
ground fitting residuals in limiting the sensitivity of our signal ex- 
traction using the skewness. For other statistics, great care may be 
required to model and remove the foregrounds to high accuracy 
( [Morales, Bowman & HewiT3l2006h . Dealing with polarized fore- 
grounds is also a concern, which will be approached in future work. 

Clearly, it is also desirable to model the CS itself accurately, 
especially when using higher order statistics. A good extraction 
scheme should work for a wide range of reionization scenarios, 
and ideally should be able to distinguish between them. We have 
therefore tested our scheme with three detailed models in which 
the ionization history is quite different. Though in each case we 
see the skewness become relatively small before rising to large val- 
ues - behaviour which may be generic ~ a more extended period 
of reionization stretches out these features. It is possible, however, 
that more exotic sources which we have neglected would cause dif- 
ferent behaviour. If, for example, decaying dark matter makes a 
significant contribution to heating the intergalactic medium before 
reionization or during its early stages, we can expect this heating to 
be uneven, the rate of energy deposition depending on the square of 
the density. Then we can no longer assume that the hydrogen spin 
temperature, Ts, is much larger than the CMB temperature, Tcmb, 
everywhere. We would have to multiply Equation[3]by a position- 
dependent factor (Ts — Tcmb) /2s, which could result in non-zero 
skewness even if the neutral fraction is approximately unity every- 
where. 

The main limitation of these simulations when it comes to 
testing our extraction scheme, however, is their size. To generate 
maps with the area of one LOFAR EoR window we must tile our 
datacubes in the image plane. In some cases this may be unreal- 
istic: when the size of individual ionized bubbles becomes com- 
parable to the size of the simulation box, a slice through the box 
can no longer be considered to be a representative slice of the Uni- 
verse. We have argued that this may not be important for one-point 
statistics (and if anything, having a larger number of independent 
volumes contributing to each image would improve our signal to 
noise ratio). Firstly, each pixel is in any case affected differently 
by foregrounds and noise which are much larger than the CS. Sec- 
ondly, nearby frequency slices are at a similar stage of reionization 
but may otherwise be sufficiently weakly correlated that smooth- 
ing along the frequency direction after extraction can help recover 
a clearer trend. Spatial statistics will clearly be more seriously af- 
fected by tiling. Note that there are also periodic repetitions in the 
frequency (redshift) direction in the simulated CS. This can be seen 
in the high redshift portion of the curve corresponding to the evo- 
lution of (STh) in the f250C simulation in Fig. |3] The onset of 
reionization appears to break this periodicity somewhat: for exam- 
ple, the mean ionized fraction can change significantly between two 
redshifts separated by a comoving radial distance corresponding to 
the size of the simulation box. It therefore seems to be no obstacle 
to robustly recovering the overall trends. 



6 SUMMARY AND CONCLUSIONS 

Many statistics have been put forward to characterize the 21cm 
emission from the EoR, the power spectrum probably being 



emission trom tne toK, the power spectrum probably being 
the most frequently studied jBarkana I2OO8I : iLidz et al.l bOOsI: 
iPritchard & Loebll2008l : ISethi & Haimaij|2008l to choose some re- 
cent examples). We have suggested that higher-order statistics may 
be useful not only to characterize a CS cube that has been cleaned 
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of foregrounds, noise and instrumental effects, but also to extract 
the signature of reionization from these corrupting influences in 
the first place. The skewness of the one-point distribution of bright- 
ness temperature, measured as a function of observed frequency (or 
equivalently as a function of redshift), is one such promising statis- 
tic. 

The three detailed simulations of reionization which we have 
studied show a strong evolution of the skewness with redshift. 
Some of the features of this evolution appear to be generic and can 
be readily understood: in the early stages of reionization the skew- 
ness drops below that of the underlying density field as the first ion- 
ized bubbles, from which the emission is negligible, are formed. As 
reionization progresses, the majority of the volume becomes ion- 
ized and the skewness increases again, becoming very large at low 
redshift when the distribution of brightness temperature is peaked 
at zero, with a tail extending to large values. In simulation f250C 
there is a well defined dip in the skewness with a width A2 « 1. 
Our other simulation in which the Universe is reionized entirely by 
stars (T-star) shows a more gradual change, with the epoch of reion- 
ization extending throughout the redshift range probed by LOFAR. 
A third simulation, T-QSO, in which QSOs reionize the Universe, 
shows an intermediate behaviour. 

By combining these simulations with models of the fore- 
grounds, noise and instrumental response, we have generated dat- 
acubes which are intended to simulate the output of the LOFAR 
EoR experiment. We have studied two cases: firstly, one in which 
we smooth the foregrounds and signal to the resolution of the tele- 
scope using a Gaussian kernel, then add uncorrelated Gaussian 
noise; secondly, one in which we degrade the foregrounds and noise 
to the resolution of the telescope using a realistic PSF, and add 
noise which is uncorrelated in the Fourier plane rather than the im- 
age plane, producing what we refer to as 'dirty' images. In the for- 
mer case, we can see the signature of reionization in the skewness 
by fitting out the foregroimds to obtain residual images, and then 
denoising these images with a simple smoothing operation. The 
skewness in these images as a function of redshift shows significant 
evidence of reionization. The result is quite robust to the details 
of the foreground fitting and the smoothing. Under- or over-fitting 
the foregrounds affects the recovered skewness less severely than 
the recovered variance. Extracting a signal from the dirty cubes 
requires a more sophisticated denoising scheme. In an optimistic 
scenario where the correlation matrices of the original signal and 
of the noise are known, we can again recover the evolution of the 
skewness quite cleanly using Wiener deconvolution. 

We have touched upon some areas for improvement: simula- 
tions which remain realistic but extend to larger scales and exhibit 
an even greater range of reionization histories; taking into account 
the polarization of the foregrounds and the instnmiental response, 
and incorporating new observational constraints as they arrive; test- 
ing the minimal assumptions we must make about the signal in or- 
der for our extraction scheme to work, for example whether a poor 
estimate of the correlation matrix of the CS seriously affects the ex- 
tracted skewness; and studying a wider range of statistics beyond 
the variance and power spectrum. All of these will be areas for 
future work. Even at this stage, however, our results justify some 
optimism that the new generation of radio telescopes can detect the 
signature of reionization using higher-order statistics. 
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